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Abstract 

The saturation conditions for bending modes in inhomogeneous thin stellar disks that 
follow from an analysis of the dispersion relation are compared with those derived from 
A^-body simulations. In the central regions of inhomogeneous disks, the reserve of disk 
strength against the growth of bending instability is smaller than that for a homogeneous 
layer. The spheroidal component (a dark halo, a bulge) is shown to have a stabilizing 
effect. The latter turns out to depend not only on the total mass of the spherical compo- 
nent, but also on the degree of mass concentration toward the center. We conclude that 
the presence of a compact (not necessarily massive) bulge in spiral galaxies may prove to 
be enough to suppress the bending perturbations that increase the disk thickness. This 
conclusion is corroborated by our A^-body simulations in which we simulated the evolu- 
tion of almost equilibrium, but unstable finite-thickness disks in the presence of spheroidal 
components. The final disk thickness at the same total mass of the spherical component 
(dark halo + bulge) has been found to be much smaller than that in the simulations where 
a concentrated bulge is present. 
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1 INTRODUCTION 



A peculiarity of the disks in spiral galaxies is that these are rather thin objects. Their 
thickness is several times smaller than the radial scale length. How far stars can go from 
the principal galactic plane due to their vertical random velocity component determines 
the disk thickness at fixed star surface density. The larger the velocity dispersion, the 
thicker the disk. The random velocities of young stars are known to be low, but the 
stellar ensemble can subsequently heat up through various relaxation processes; i.e., the 
random velocity dispersion can increase. Thus, the stellar disk thickness depends on how 
effective the relaxation processes are in galaxies, and it is ultimately determined by the 
factors that suppress or trigger the various heating mechanisms. Three basic stellar disk 
heating mechanisms are commonly discussed in the literature: the scattering of stars by 
giant molecular clouds (Spitzer and Schwarzschild 1951, 1953), the scattering by transient 
spiral density waves (among the first results of numerical simulations are those obtained 
by Sellwood and Carlberg (1984)), and the heating of the ensemble of stars that constitute 
the disks of spiral galaxies as they interact with external sources, for example, with low- 
mass satellites (see, e.g., Walker et al. 1999; Velasquez and White 1999). 

A second remarkable peculiarity of the stellar disks is that their structure is unusually 
"fragile" . This peculiarity was revealed by a linear analysis of the coUisionless Boltzmann 
equation and has been repeatedly illustrated by A'"-body simulations. Numerous studies 
have shown that the initially regular structure of the stellar disks can change radically 
due to the growth of various instabilities, which give rise to large-scale structures both 
in the disk plane (bars, spiral arms, rings) and in the vertical direction (warps). The 
analytically and numerically obtained saturation conditions for unstable modes impose 
the most severe restrictions on the global structural and dynamical parameters of the 
stellar disks. 

A local analysis suggests that the stars at a given distance R from the disk center must 
have a radial velocity dispersion (Tr(-R) larger than some minimum critical value o'^(-R) 
(Toomre 1964) for the disk to be gravitationally stable in the region under consideration 
(at least against the growth of axisymmetric perturbations). In addition, for the disk to 
be stable against the growth of bending perturbations, the ratio of the vertical and radial 
velocity dispersions Ozjo^ must also be larger than some value approximately equal to 
0.2 — 0.37 (Toomre 1966; Kulsrud et al. 1971; Polyachenko and Shukhman 1977; Araki 
1985). The latter quantity determines the minimum thickness of the galactic disk, and 
its comparison with the observed value allows the contribution of the bending instability 
to the vertical disk heating to be estimated and compared with the contribution of other 
relaxation mechanisms. 

On the other hand, if we exclude the heating mechanisms mentioned above from our 
analysis and take, as is commonly done, the condition for marginal disk stability against 
the growth of bending modes by fixing u^/ur at a level of the linear approximation, 
0.29 — 0.37, then the velocity dispersion (T^ at the same star surface density will decrease 
with decreasing aR (see, e.g., Zasov et al. 1991). In this case, the disk will have a smaller 
thickness. 

The presence of a spheroidal component, for example, a dark halo is known to produce 
a stabilizing effect and to decrease the minimum value of cr^ required for gravitational 
stability. Consequently, the disks embedded in a massive halo, on average, must have low 
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values of cr^ and be, on average, thinner. Based on similar reasoning, Zasov et al. (1991, 
2002) showed that the relative disk thickness z^/h {zq is the half-thickness of the disk, 
and h is the radial exponential scale length) is proportional to M^/Mx,, where and 
Mt are, respectively, the mass of the disk and the total mass of the galaxy within a fixed 
radius. These authors also concluded that a small disk thickness suggests the existence 
of a massive dark halo in the galaxy. Moreover, based on A^-body simulations, Zasov et 
al. (1991) and Mikhailova et al. (2001) constructed a dependence that allows the relative 
mass of the dark halo M\^/Ma to be estimated from Zq/H. However, as we show below, the 
relationship between the relative disk thickness and the mass of the spheroidal component 
is more complex. 

In this paper, we numerically analyze the saturation conditions for bending instability 
in inhomogeneous three-dimensional stellar disks at nonlinear stages in the presence of 
a spheroidal component of different nature (a stellar bulge and a dark halo) and the 
constraints imposed on the final disk thickness. 

2 PECULIARITIES OF THE GROWTH OF 
BENDING INSTABILITY IN INHOMOGE- 
NEOUS THIN STELLAR DISKS 

2.1 Global Modes 

To understand how the growth of bending instability in inhomogeneous disks differs from 
that in homogeneous disks, let us first turn to the result of Toomre (1966). Toomre was 
the first to derive the dispersion relation for long-wavelength bending perturbations in an 
infinitely thin gravitating layer with a nonzero velocity dispersion 

uo"^ = 2TiGT.\k\ - alk'^ . (1) 

where S is the star surface density of the layer, and is the velocity dispersion along 
a particular coordinate axis in the plane of the layer. It follows from Eq. that the 
perturbations with a wavelength A = 27r//c > Aj = a^/GS are stable, since tu^ > in this 
range. 

Relation was derived for an infinitely thin disk and, when applied to finite-thickness 
disks, is valid only for perturbations with a wavelength longer than the vertical scale height 
of the system, i.e., for A >> zq, where Zq is the half-thickness of the layer. It can be shown 
from general considerations that the longest-wavelength perturbations in finite-thickness 

(Jx 

disks with a wavelength \ < \2 ^ Zq — are stable, since they are smeared by thermal 

motions in the plane of the layer (see, e.g., Polyachenko and Shukhman 1997). Having 
derived the corresponding dispersion relation, Polyachenko and Shukhman (1977) were 
the first to find the exact location of the stability boundary in the short-wavelength range 
for a homogeneous fiat finite-thickness layer. Araki (1985) (see also Merritt and Sellwood 
1994) obtained a similar result for a homogeneous layer with a vertical density pro?le 
close to the observed one in real galaxies^ p{R, z) = p(i?, 0)sech^(z/2;o). As regards the 

^This profile corresponds to the model of an isothermal layer (Spitzer 1942) and describes well the 
vertical density variations observed in galaxies (van der Kruit and Searle 1981). 
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intermediate-wavelength (A2 < A < Aj) perturbations, they are unstable, as follows from 
the results of the studies mentioned above. 

As the disk thickness Zq increases, the wavelength A2 increases and tends to Aj. When 
-^2 = A J, the disk is stabilized against bending perturbations of any wavelengths. The 
following analytical estimate in the linear approximation obtained both from qualitative 
considerations (Toomre 1966; Kulsrud et al. 1971) and from an accurate analysis of the 
dispersion relation for a finite-thickness layer (Polyachenko and Shukhman 1977; Araki 
1985) is valid: 

((T,/a,),,^ 0.29 - 0.37. (2) 
The instability is completely suppressed if az/dx > {cTz/cTx)^^ and grows if Ozjox < 

As regards the inhomogeneous models, the first thing that radically distinguishes the 
growth of bending instability in inhomogeneous disks from that in a homogeneous layer 
is the existence of global unstable bending modes with a wavelength longer than the disk 
scale length. This conclusion follows from an analysis of the equation that describes the 
evolution of long-wavelength bending perturbations in an in infinitely thin disk with a ra- 
dially decreasing density (Polyachenko and Shukhman 1977; Merritt and Sellwood 1994). 
For special disk models, it was shown that the region of stable long-wavelength pertur- 
bations narrows significantly in this case (see, e.g.. Fig. 2 from Merritt and Sellwood 
(1994)). This is attributed to the fact that the restoring force from the perturbation that 
grows in an inhomogeneous disk (Merritt and Sellwood 1994) or in a radially bounded 
disk (Polyachenko and Shukhman 1977) is always weaker than the corresponding restor- 
ing force in a homogeneous infinite layer. This fact was demonstrated more clearly by 
Sellwood (1996), who noted that the dispersion relation ((T]) could be used to analyze the 
bending instability in inhomogeneous disks (at least qualitatively) if another term related 
to the restoring force from the unperturbed disk is added to it: 

uj^ = ul + 27fG (3) 

where z/j = ^yd'^^^iR, z)/dz'^ is the vertical oscillation frequency of the stars, and $d(-R, z) 
is the potential of the disk. 

For disks with a nonfiat rotation curve, the additional term z/^ can play a destabilizing 
role. As was noted by Sellwood (1996), for an infinitely thin inhomogeneous disk. 



52$d(i?, z) 



~ dz^ 
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' (4) 



z=0 ^ 



(fc,d is the circular rotational velocity of the disk), and < for dv^^^/dR > 0. 

Thus, an additional expulsive force from the unperturbed disk emerges in the central 
regions where the rotation curve rises. The destabilizing effect of the disk in the central 
regions gives rise to another region at small wave numbers (long A) where cu^ < 0. This 
region is responsible for the growth of largescale bending instability and the emergence 
of global modes. In this case, one might expect a larger disk thickness and a larger 
value ctz/ctx- (o"z/o"r) than those in homogeneous models to be required to suppress the 
instability. 

Applying all of this reasoning to finite-thickness disks is not quite obvious. On the 
one hand, it may be assumed that the restoring force from the unperturbed disk Fz = 
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—d^d/dz for i? — i> behaves as ~ —GM(iz/\z\^ (here, is the total mass of the 
disk) starting from some z, i.e., decreases (in magnitude) with increasing z. Consequently, 
z/J = d'^^d/dz'^ = —dFz/dz becomes negative (Sellwood 1996). 

On the other hand, we can exactly calculate the z dependence of z/^ for a given R from 
the general formula 

2^ \ ^T? m ^ f Gpd{r'){z' - z) 3 , 
^d(r) = -dFjdz = -Q^J ^737p ' 

for the density profile that is commonly used to describe the disks of spiral galaxies: 

where h is the exponential scale length of the disk, zq is the vertical scale height, and 
Md is the total mass of the disk. The derived dependence^ for various R is shown in 
Fig. Figure CJi shows the z dependence of the vertical force Fz, which necessarily has 
an extremum at some z. In the central regions, is negative at z^2zo (on the periphery, 
the passage to the negative region occurs at even larger z/zq). This implies that all of 
the above reasoning for finitethickness disks is directly applicable only to largeamplitude 
perturbations. It is also clear that the bellshaped (or axisymmetric) mode with the 
azimuthal number m = must be most unstable in the central regions. If the amplitude 
of the axisymmetric bend increases significantly during the growth of bending instability, 
then it can raise the central stars above the plane of the disk and bring them into a 
"dangerous" zone where v\ < 0; subsequently, the self-gravity of the disk itself will 
contribute to the growth of a bend to the point of saturation. 



2.2 The Central Regions of Hot Infinitely Thin Disks 

The second peculiarity of the growth of bending instability in inhomogeneous infinitely 
thin disks is as follows. The degree of disk instability against bending perturbations 
depends on the degree of disk "heating" in the plane. Exact and approximate analyses 
of the dispersion relation for an inhomogeneous infinitely thin disk (Polyachenko and 
Shukhman 1977; Merritt and Sellwood 1994) show that the modes with increasingly 
large azimuthal numbers in the central regions become unstable for all wavelengths as 
the fraction of the kinetic energy contained in the random motions in the disk plane 
increases. The conclusion about the existence of a Toomre parameter Qt (Toomre 1964), 
which characterizes the degree of heating of the stellar disk in the plane, at which the 
disk cannot be stabilized against bending perturbations of any wavelengths can also be 
drawn from Eq. i.e., the equation written for long- wavelength perturbations in a 
homogeneous layer, but "rectifided" by the term i/J to include the inhomogeneity effects. 

^We see from Eq. (O that triple integrals over infinite intervals must generally be calculated to 
determine v^. This can be easily done by using adaptive algorithms for calculating the integrals. In 
our integration, we used the gsl library (information about the gsl (GNU Scientific Library) project can 
be found at http://sourses/redhat/com/gsl). For models with known analytical density-potential pairs, 
for example, the Miyamoto-Nagai disk (see Binney and Tremaine (1987), p. 44), a comparison of the 
analytical z dependence of v'^ at a given R with that calculated numerically by triple integration yielded 
a close match. 
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Indeed, u"^ < for any k if 



2 (ttGE)^ 



,,2 _ _ ^2 
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Given (jH), 



where f2d = is the angular velocity, and = ^^-y I 1 H --^ is the square of 

the epicyclic frequency. Then, 

u'<Q 2fi^-«:2+ (![^<0. (7) 

For the central regions of a rigidly rotating disk {n\ = 4^2^), we obtain from ((7|) 

or 

cu^ < Qt > "\/2 for any wavelength, (8) 

where Qt = <^x/<^'r is the Toomre parameter^ (Toomre 1964). 

Condition (jHl) is not an exact criterion; it is only an estimation relation. However, 
it shows that the central regions of hot stellar disks (Qt >> 1) with a large reserve of 
strength against the growth of instabilities in the plane of the disk (bars, spiral arms) 
cannot be stabilized against the growth of bending perturbations of any wavelengths. 
The theory constructed for a homogeneous infinitely thin layer does not yield this regime. 
It should be borne in mind, however, that the contribution of the destabilizing term (z/J) 
must decrease with increasing disk thickness; therefore, the instability will be saturated 
at a large, but finite disk thickness. However, the following might be expected: other 
things being equal, the hotter the initial model in the plane, i.e., the larger the Toomre 
parameter Qt, the higher the saturation level. 



3 BENDING INSTABILITY: NUMERICAL SIMU- 
LATIONS OF THREE-DIMENSIONAL DISKS 

The nonlinear growth stages of bending instability in inhomogeneous finite-thickness stel- 
lar disks have been extensively investigated by numerically solving the gravitational A^- 
body problem for various stellar disk models. Raha et al. (1991) first observed the bending 
instability of bars in their numerical simulations; Sellwood and Merritt (1996) and Merritt 
and Sellwood (1994) studied the nonlinear regime of bending instability in nonrotating 
disks with the radial density profiles that corresponded to the Kuzmin-Toomre model 
(see, e.g., Binney and Tremaine (1987), p. 43); Griv et al. (1998) numerically analyzed 
the development of a bend in a layer of newly formed stars; Tseng (2000) simulated the 
evolution of the vertical structure of a homogeneous, initially thin disk of finite radius; 

■^Actually, the Toomre parameter should have been defined as ctt^/ct^', but the difference does not 
matter in our case. 
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Sotnikova and Rodionov (2003) considered a rotating disk with an exponential density 
profile along the R axis and assumed the presence of a dark halo in the system. The 
latter authors analyzed the question of how the evolution of initially equilibrium thin 
disks depends on the governing parameters of the bending instability, which include the 
initial disk half-thickness zq, the Toomre parameter Qt, and the relative mass of the 
dark halo M^/M^ within a fixed radius. Below, we hst the most important conclusions 
that follow from the A^-body simulations described in the literature. These conclusions in 
many respects agree with those given in the section entitled "Peculiarities of the growth 
of bending instability. . . " of this paper for inhomogeneous infinitely thin disks. 

(1) In inhomogeneous models, all of the experimentally observed modes are global, 
i.e., the scale length of the unstable perturbations is larger than the disk scale. The linear 
theory constructed for homogeneous models yields no such result. 

(2) The saturation level for the bending instability depends on cr^ (or Qt)- The 
larger the reserve of disk strength against perturbations in the disk plane, the greater 
the dificulty to stabilize the disk against the growth of bending perturbations. A rapid 
and significant (occasionally catastrophic) increase in the disk thickness, particularly in 
the central regions, to values that are severalfold larger than those yielded by a linear 
analysis for homogeneous, moderately hot models (Toomre 1966; Kulsrud et al. 1971; 
Polyachenko and Shukhman 1977; Araki 1985) was observed in all of the simulations 
with initially hot disks (Sellwood and Merritt 1994; Tseng 2000; Sotnikova and Rodionov 
2003). This mechanism may be responsible for the formation of central bulges in spiral 
galaxies, at least it can feed the spherical component with new objects. 

(3) The central regions of the disk are most unstable (Sellwood and Merritt 1994; 
Merritt and Sellwood 1994; Griv and Chiueh 1998; Griv et al. 2002; Sotnikova and Rodi- 
onov 2003). It is here that the bending modes are formed. Subsequently, their amplitude 
increases, sometimes significantly. This is particularly true for the perturbations with 
m = 0. The nonaxisymmetric modes (with the azimuthal numbers m = 1 and m = 2) 
drift to the disk periphery, temporarily creating the effect of a largescale warp of the 
entire galaxy, and are then damped''. The instability in the central regions of the stellar 
disks is saturated at (Sotnikova and Rodionov 2003) 

cTz/ctr ^ 0.75 - 0.8 . 

(4) The presence of a massive dark halo, which was included in the models by Sotnikova 
and Rodionov (2003), has always been a stabihzing factor that suppresses the growth of 
bending modes. This effect appears to have been first described qualitatively by Zasov et 
al. (1991). 

^In contrast to the result of Griv et al. (2002), the largescale S-shaped or U-shaped warp of the 
galactic edge always disappeared on long time scales (l 5 Gyr) in all of the simulations by Sotnikova and 
Rodionov (2003). 
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4 THE STABILIZING EFFECT OF THE 
SPHEROIDAL COMPONENT: A QUALITA- 
TIVE ANALYSIS 

The remarkable agreement (at least on a qualitative level) between the conclusions that 
follow from the analysis of the dispersion relation for an infinitely thin disk and the 
results of numerical simulations with three-dimensional disks allows us to analyze the 
applicability of yet another conclusion that can be drawn from Eq. (jH)) to finite-thickness 
disks. Since the central regions of the disk are most unstable, it is important to separate 
out the factors that have a stabilizing effect precisely on these regions. An additional 
spherical component (a dark halo, a bulge) can be such a stabilizing factor. In this case, 
another term related to the restoring force exerted from the spherical component appears 
in Eq. with 



sph 



>0. (9) 



This term was also introduced by Sellwood (1996) when analyzing the dispersion 
relation for the longwavelength bending perturbations of an infinitely thin disk, but its 
role was not studied. 

Zasov et al. (1991, 2002) and Mikhailova et al. (2001) concluded that the minimum 
possible relative thickness of an equilibrium stellar disk, z^/h, decreases with increasing 
relative mass of the dark halo. Let us consider the relationship between the stellar disk 
thickness and the mass of the spheroidal component in terms of the stabilization conditions 
for the bending modes in inhomogeneous thin disks. 

We take specific bulge and halo models (these models were subsequently used in our 
numerical calculations). A Plummer sphere is taken as the bulge model. Its potential is 
(see, e.g., Binney and Tremaine (1987), pp. 42-43) 

(r^ + ag) ' 

where Mb is the total mass of the bulge, and ab is the scale length of the matter distri- 
bution. We describe the potential of the dark halo in terms of the logarithmic potential 
(see, e.g., Binney and Tremaine (1987), p. 46) 

*h(r) = ^ln(r2 + a^) + const, (11) 

where Oh is the scale length, and foo is the velocity of a particle in a circular orbit of 
infinite radius. The parameter f oo is related to the mass of the halo with a sphere of given 

radius r by Mh(r) = — k). 

G + 

For the additional stabilizing term in the dispersion relation, the models of the spher- 
ical components (the bulge and the halo) that we used yield 



(i?2 + a2)3/2 ' ^ R^ + al 
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The stabilizing effect of the spherical component weakens at large R, but the disk itself 
in the peripheral regions has a stabilizing effect: ~ —GM^z/R^ at \z\ « R ^ = 
—dFz/dz > 0. On the other hand, the strength of the effect increases in the central (most 
unstable) regions (i.e., for R —>■ 0); this strength depends not only on the total mass of the 
spherical component (Mb or foo), but also on the degree of matter concentration toward 
the center, Ob and ah. It thus follows that the presence of a compact (not necessarily 
massive) bulge in galaxies may prove to be enough to suppress the bending perturbations. 
This implies that the disks of galaxies with compact bulges can be as thin as the disks 
embedded in a massive dark halo. 

To test our conclusion, we carried out a series of numerical simulations. 

5 THE STABILIZATION OF BENDING PERTUR- 
BATIONS BY A COMPACT BULGE: iV-BODY 
SIMULATIONS 

5.1 The Method 

We used an algorithm based on the hierarchical tree construction method (Barnes and Hut 
1986) to simulate the evolution of a self-gravitating stellar disk. In our calculations, we 
always included the quadrupole term in the Laplace expansion of the potential produced 
by groups of distant bodies. The parameter 9 (Hernquist 1993) that is responsible for the 
accuracy of calculating the gravitational force was chosen to be 0.7 in all our simulations. 
The NEMO software package (http:/ /astro. udm.edu/nemo; Teuben 1995) was taken as 
the basis. We enhanced the capabilities of this package by including several original codes 
for specifying the equilibrium initial conditions in a fiat stellar system^ and supplemented 
it with new codes that allow us to easily analyze the data obtained and to present them 
in convenient graphical and video formats. 

5.2 The Numerical Model 

When constructing the galaxy model, we distinguished two components in it: a self- 
gravitating stellar disk and a spherically symmetric component that was described in 
terms of the external static potential, which is a superposition of two potentials, ()10|) for 
the bulge and (jllj) for the dark halo. At large distances R from the center of the stellar 
system, in the region where the halo dominates, potential (fTTj) yields a fiat rotation curve. 
The disk was represented by a system of gravitating bodies with the density profile (jHl). 

The initial conditions in the A^-body problem suggest specifying the mass, position 
in space, and three velocity components for each particle. The particle coordinates are 
naturally determined in accordance with the disk matter density profile the distant 
regions of the disk are disregarded. We took only those particles for which the cylindrical 
radius R < -Rmax and 1^1 < -Zmax- The mass of all particles was assumed to be the 

^The mkexphot code for specifying equilibrium stellar disk models from the NEMO package has 
limitations on the parameters of the outer halo and does not enable the gravitational field of the bulge 
to be specified. The models constructed using our codes closely agree with those obtained using the 
mkexphot code for identical parameters. 
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same. The total mass of the particles was equal to the mass of the disk region under 
consideration (i.e., the disk region for which R < -Rmax and \z\ < -Zmax)- The particle 
velocities were specified using the equilibrium Jeans equations by the standard technique 
(see, e.g., Hernquist (1993), Section 2.2.3). 



5.3 Specifying the Velocity Field in the Model Galaxy 

To specify the initial particle velocities for a disk that is in equilibrium in the plane and 
in the vertical direction, we make the following assumptions: 

(1) The velocity distribution function is the Schwarzschild one; in other words, the 
particle velocity distribution function has only four nonzero moments: the mean azimuthal 
velocity v^p, the radial velocity dispersion 0"^?, the azimuthal velocity dispersion cr^, and 
the vertical velocity dispersion^ a^- 

(2) All four moments depend only on the cylindrical radius R and do not depend on 

z. 

(3) The epicyclic approximation is valid^. 

(4) cr^ is proportional to the surface density of the stellar disk, i.e., oc exp (—R/2h) 
(this assumption agrees well with the observational data; see, e.g., van der Kruit and 
Searle 1981). 

The following relations for the moments (in which the R dependence was omitted for 
simplicity) can then be derived from the Jeans equations (see, e.g., Binney and Tremaine 
1987): 

(12) 

where Vc is the circular velocity of a particle placed in the total potential of the disk and 

the spherical component, the bulge and the halo, (f^ = "J^cd + ^cb + ^ch! "^cb — ^^w'l 

' [ '_ ' oR 

.d^h. „ Vr . , , , . , /„ t>2 1 dv'i 



V = R——), = — is the angular velocity, and n = \ 2-^ + — —p- is the epicyclic 

oR R \ R^ R dR 

frequency. 

The circular velocities for the bulge ()1U|) and the halo (fTTj) have analytical expressions. 
The circular velocity for the disk (jH)) can be determined by numerical integration (see the 
section entitled "Global modes" ) using the general formula 



fr' - r) ■ R 
r| 



where R is the projection of the vector r onto the disk plane. 

The Jeans equations that are used to derive relations (|T2|l are known to provide no 
exact disk equilibrium (see, e.g., Binney and Tremaine 1987). Moreover, the last relation 
in ()12j) . which follows from the vertical equilibrium condition for a disk with the density 

^As many members of the astronomical society, we have the bad habit of calhng the standard of the 
distribution function the dispersion. 

^In the central regions of the disk, this approximation breaksdown. 
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profile (jni) and a independent cr^, was written without including the influence of the 
additional spheroidal components. The adjustment to equilibrium occurs on time scales 
of the order of several vertical oscillation times This time was always no longer than 

100 — 120 integration time steps for the equations of motion (the thinner the disk, the 
shorter this time) and much shorter than the instability growth time scale in the disk. 
Moreover, in the context of the problem of the growth and saturation of unstable modes, 
a small deviation of the disk from equilibrium at the initial time may be treated as an 
additional initial perturbation. 

The fourth assumption (see above) about the radial velocity dispersion cr\{R) causes 
diflculties in calculating in the central regions. In the flrst equation of system [T21 v"^ 
is occasionally negative at small R (since cr^(-R) rapidly increases toward the center, the 
last term on the right-hand side can make a large negative contribution). For this reason, 
the dependence for was reduced at the center (Hernquist 1993): 

an oc exp [-^R' -2al/2h) . (14) 

If the parameter is taken to be h/A — h/2, then this proves to be enough to properly 
calculate v^^. In the central regions of the disk, adjusted to in such a way that 
the ratio OzjoB. was constant at a given half-thickness = const in the initial model 
throughout the disk. 

The proportionality factor in can be determined via the Toomre parameter Qt 
at some radius R^-ei- 

aniRrei) = Qt al(i?ref) = Q^ 3.36- Sd(/g.ef) _ ^^^^ 

The sought proportionality factor can be obtained from ((T3j) and |TKj) . Specifying 
Qt; which ensures disk stability in the plane, at i^^ef ~ 2.5h yields the condition 
Qt{R) > QT(-Rref) (Hernquist 1993). The latter, in turn, ensures a stability level against 
perturbations in the disk plane no lower than that at -Rref- The initial half-thickness zq 
for the adopted QT(-Rref) was chosen in such a way that ?z/?R was less than 0.3 — 0.4, 
which ensured initial instability against the growth of bending modes. 



5.4 Parameters of the Problem 

All of the results discussed below are presented in the following system of units: the 
gravitational constant is G = 1, the unit of length is = 1 kpc, and the unit of time 
is t„ = 1 Myr. The unit of mass is then M„ = Rl/Gtl = 22.2 ■ IO^^Mq, and the unit of 
velocity is f„ = Ru/tu = 978 km s~^. 

The number of bodies in the simulations was = 300 000 (in several cases, 500 000 
and 600 000). The force of interaction between two particles with coordinates rj and 
and masses rrii and rrij was modifled, as is commonly done, as follows: 

where e is the softening length of the potential produced by an individual particle. When 
collisionless systems are simulated, this parameter is introduced for two reasons. First, the 
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divergence of the interaction force in close particle-particle encounters must be avoided 
when integrating the equations of motion. Second, when the phase density of a collisionless 
system is represented by a finite number of particles, the inevitable fluctuations in the 
particle distribution must be smoothed in such a way that the forces acting in the system 
being simulated were are to the forces acting between the particles in a system with a 
smoother density profile. The softening length e was chosen to be 0.02. This value is 
approximately a factor of 2 or 3 smaller than the mean separation between the particles 
(at = 300 000) within the region containing half of the disk mass. On the one hand, 
it matches the criterion for choosing e based on minimization of the mean irregular force 
(Merritt 1996) and, on the other hand, allows the vertical structure of thin disks to be 
adequately resolved. 

To integrate the equations of motion for particles in the self-consistent potential of 
the disk and the external field produced by the spheroidal component, we used a leapfrog 
scheme that ensured the second order of accuracy in time step. The time step was 0.5 (in 
several models, 0.25)^. 

We constructed a total of about 60 models. The entire set of models can be arbitrarily 
divided into two classes: the models with and without bulges. The scale length of the 
density distribution in the bulge ay, was assumed to be equal to 0.5 almost for all of 
the models with bulges. In several models without bulges, we chose a concentrated halo 
(tth = 2). In all of the remaining cases, the halo was "looser" (oh = 10). The total relative 
mass of the spheroidal components ^ = Msph(4/i)/Md(4/i) was varied over the range 0.25 
to 3.5 

The disk in our models has the following parameters^: h = 3.5 and the disk mass (in 
dimensional units) M^{4h) = (4 — 8) x lO^^M©. The initial thickness was varied over 
the range zq — 0.1 — 0.5. In order not to abruptly cut off the model disk at the radius 
corresponding to the optical radius (~ Ah), we chose i?max = 25, with Zmax = 5. The 
smoothing parameter of the initial radial profile of the velocity dispersion (j/j is = 1. 
The parameter in the discussion of our simulations is given for the radius i?ref — 8.5. 

5.5 Simulation Results 

Previously (Sotnikova and Rodionov 2003), we showed that there are two distinct vertical 
stellar disk relaxation mechanisms related to bending instability: the bending instability 
of the entire disk and the bending instability of the bar forming in the disk. The for- 
mer mechanism dominated in galaxies that are hot in the plane ((5t~2.0), and the bar 
formation was suppressed in this case; the latter mechanism dominated in galaxies with 
a moderate Toomre parameter Qt (such galaxies were unstable against the growth of a 
bar mode). In the simulations whose results are presented and analyzed below, we also 
considered two distinct cases: hot disks (Qt = 2.0) in which only bending instability 
developed, and cooler models (Qt = 1-5) — here, we observed the combined effect of the 
two types of instability. 

Hot disks. For hot (in the plane) disks (Qt = 2.0), we revealed distinct patterns of 
growth and saturation of bending perturbations that are consistent with the conclusions 

^Thc choice of the time step is limited above by e — the particle must take at least one step on the 
smoothing length. 

^ These parameters are close those of the disk in our Galaxy. 
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following from a qualitative analysis of the dispersion relation Q. 

The stabilizing effect from the presence of a massive spherical component that was 
discussed in the section entitled "The stabilizing effect of the spheroidal component. . . " 
is clearly seen in Fig. |2l This figure shows the radial profile of (Jz/(^r for model 12 with 
Mb = and /i = 3.0. Throughout the disk, az/cjR was set at ~ 0.35 and was determined 
not by the bending instability, but by the disk heating through the scattering of stars 
by inhomogeneities related to the different disk thickness in the spiral arms and in the 
interarm space (Sotnikova and Rodionov 2003). 

We will demonstrate the stabilizing effect of a compact (not necessarily massive) bulge 
comparable to the effect of a massive dark halo with a broader density profile using the 
results obtained for the following group of models as an example: 50, 76, 75, 49, and 53. 
In all five models, the total mass of the spheroidal component is the same and accounts 
for half of the diskmass within four exponential disk scale lengths, /i = 0.5, but it is 
differently distributed between the two spherical subcomponents. In model 50, all of the 
mass is contained in the halo (/i = /ih = Mh(4/i)/Md(4/i)); in model 53, only a compact 
bulge is present (yU = = Mb(4/i)/Md(4/i)). The remaining models are intermediate 
between the two extreme models. The initial thickness for all of the models was chosen 
to be the same, Zq = 0.1. The rotation curves for these models are shown in Fig. El The 
variety of the shown curves to some extent reflects the actual variety of rotation curves 
for spiral galaxies. 

We traced the evolution of these models up to t = 5000. Figure IH illustrates the 
variations in the dynamical parameters of the disk and cr^ the radial and vertical 
velocity dispersions calculated at R = 2h. All of the models demonstrate an initial 
increase in and a decrease in a^. Subsequently (after t ^ 1000), the latter parameter 
reaches an approximately constant value, while az for some of the models (this is primarily 
true for the model with a massive bulge) continues to slowly increase. The number of 
particles in our models and the softening length were chosen in such a way that the two- 
body relaxation time was much longer than the time scale on which we considered the 
evolution of our numerical models. The absence of heating related to numerical relaxation 
is confirmed by the behavior of and the preservation of the pattern of evolution of the 
system as the number of particles increases to = 600000. The continuing small secular 
increase in the vertical velocity dispersion probably reflects the fact that some of our 
models did not reach a steady state^^. 

Figure El shows five frames that correspond to the late evolutionary stages of our 
model disks. As expected, the saturation level for the bending instability in model 50 was 
very high and did not match the standard linear criterion. The galaxy greatly thickened 
at the final evolutionary stages. However, when we transferred 50% of the mass from the 
halo to the bulge (model 49: yUh = yUb = 0.25), the picture changed. The saturation level 
for the bending instability became much lower. At the final evolutionary stages, the disk 
was much thinner than that in model 50. In model 53, when we placed all of the mass of 
the spherical component in the bulge, the amplitude of the observed bend was very low, 
and the galaxy remained quite thin even at the late evolutionary stages. 

^"in model 53, a bulge with a mass equal to half of the disk mass and a scale length of 500 pc is atypical 
of real galaxies. We consider this as a limiting case. 

^^If the system has no third integral of motion, then its evolution to equilibrium must eventually lead 
to the relation cr^ — an. 
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The disk thickness can be quantitatively estimated as the root-mean-square (rms) 
value of the z coordinates of the disk particles, Zyms = \/ < z^ > — < z >2. This estimate 
is commonly encountered in the literature. It can be shown that for the vertical density 
profile (0), the relationship between this parameter and zq is given by Zrms = 'it/2\/3zo ^ 
0.9l2o- In practice, however, this parameter proved to be a not very good characteristic 
of the thickness. First, the fluctuations in this parameter along R were found to be 
great even when using a large number of particles if only no averaging is performed in 
concentric rings of large width. Second, the thickness calculated in this way turns out to 
be systematically overestimated due to the existing of a significant tail of the particles 
that went far from the disk plane. For these reasons, we estimated the disk thickness at 
a given distance R through the median of the absolute value of z that was designated as 
2:1/2. Twice the value of 2:1/2 is nothing but the disk thickness within which half of the 
particles is contained. For the density profile 22:1/2 = 2:0 ■ ln3 ~ l.l2:o. 

Figure ini shows the differences between the radial disk thickness profiles for model 53 
obtained by the two described methods (the averaging was performed in concentric rings; 
the ring width was AR = 0.4). We see that ^1/2 behaves much more smoothly (we have 
in mind the overall monotonic dependence of the density and the fluctuation level) than 
does the rms value of the z coordinate commonly used to estimate the disk thickness in 
A^-body simulations. 

Figure m shows the radial disk thickness profile for models 50, 76, 75, 49, and 53 at the 
time t = 3000. Note that the thickness profile for model 50 is rather unusual in shape. 
This shape is most likely attributable to the existence of X-shaped stationary orbits in 
the central regions that arise at a certain disk thickness when there are conditions for 
the resonance between the stellar oscillation frequencies in the disk plane and in the 
vertical direction (see, e.g., Patsis et al. 2002). We see the following from Fig. as 
well as from Fig. El which show the edge-on views of the model galaxies: the thinner the 
galactic disk, the larger the mass of the spheroidal component contained in a compact 
bulge. Since not all of our models reached a steady state, their thickness continues to 
slowly increase (Fig. ISJ, but the differences in thickness are always preserved. Similar 
results were obtained for all of the remaining such models with the same mass of the 
spherical component in the range Msph(4/i) = 0.25Md(4/i) to Msph(4/i) = 3.5Md(4/i). The 
stabilizing effect of a bulge was particularly pronounced in those cases where the bulk of 
the galactic mass was contained in the disk. 

The final disk thickness at fixed initial Qt "was determined only by the relative mass 
of the spheroidal component and the contribution of the bulge to this mass and did not 
depend on how far from stability the initial state of the disk was chosen. The start from 
different initial disk thicknesses led to models without any systematic differences between 
them. Figure ini illustrates this result, which is similar to that obtained in their numerical 
simulations by Sellwood and Merritt (1994). 

Thus, our three-dimensional calculations are in good agreement with the conclusion 
following from our analysis of the dispersion relation for a thin disk that a bulge is an 
effective stabilizing factor during the growth of bending instability. Moreover, since the 
initial bend is formed in the most unstable central part of the galaxy (Sotnikova and 
Rodionov 2003), it is the central regions that must be stabilized. This does not require a 
massive dark halo; the presence of a compact spherical component like a bulge will suffice. 
This suggests that the final thickness of the model galaxy depends not only on the total 
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mass of the spherical component, but also on the mass distribution in it. 

Barred galaxies. The bending instability of bars is an effective vertical disk heat- 
ing mechanism for galaxies unstable against the formation of a bar (Raha et al. 1991; 
Sotnikova and Rodionov 2003). In contrast to the warp in the entire disk, the warp 
in the bar is formed not in the central regions, but in the entire bar simultaneously 
(this is seen particularly clearly in our color two-dimensional histograms of the warp ac- 
cessible at http: / /www.astro. spbu.ru/ staff /seger/articles/warps_2002 / fig6_web.html and 
[http: / /www.astro. spbu.ru/staff/seger / articles / warps_2002 / fig7_web. html ) . Therefore, 
the conclusion that a compact bulge during the growth of bending instability in a bar 
will have the same effective stabilizing effect as that for hot stellar disks is not obvious in 
advance. 

In our simulations with Qx = 1-5 and /i^Sl-O, the warp in the bar was formed early, 
at t ~ 800. The presence of a compact bulge eventually led to the formation of thinner 
disks, although the effect itself was fairly complex. 

Within R < 1.5h, the most prominent features of the bars at late evolutionary stages 
were X-shaped structures. If the disk is viewed edge-on, then they manifest themselves 
as a bulge with an appreciable extent in the z direction with boxy isophotes. In the 
region R < 1.5h where the bar dominated, we failed to reveal any distinct patterns in 
the model disk thickness variations with increasing contribution of the bulge to the total 
mass of the spheroidal component. In general, the presence of a bulge "pushed forward" 
the bar formation time and caused the saturation level for the bending instability of a bar 
to lower. However, a further analysis is required to completely understand the processes 
during the interaction between a compact bulge and a bar. 

As regards the peripheral regions of the disks {R > 1.5h), they were always apprecia- 
bly thinner in the models with a compact bulge at late evolutionary stages. Figure ITOl 
demonstrates this effect for two groups of models with different bulge contributions: mod- 
els 52 and 51 with fi = 0.5 and models 43 and 42 with fi = 0.875. The differences in 
thickness show up most clearly in the models with a small relative mass of the spheroidal 
component. When n increases, the disk becomes very thin, as might be expected, and 
its thickness ceases to depend on how the mass is distributed between the halo and the 
bulge. 



6 CONCLUSIONS 

A comparison of the conclusions that follow from a linear analysis with the results of 
numerical simulations for three-dimensional disks shows that, in contrast to homogeneous 
models, global bending modes with a wavelength longer than the disk scale length can arise 
in inhomogeneous disks. If the amplitude of the waves during the growth of instability 
increases significantly, then they heat it up significantly in the vertical direction as they 
pass through the entire disk. Hot disks are most unstable against the growth of bending 
perturbations. An additional spheroidal component, for example, a dark halo is a factor 
that stabilizes the bending perturbations. 

Our additional qualitative analysis of the dispersion relation for inhomogeneous models 
led us to new conclusions regarding the stabilization conditions for the bending modes in 
stellar disks. These conclusions were confirmed in our numerical simulations. 
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(1) Since the central regions of the disk (particularly if the disk is hot) are most unsta- 
ble, the conditions under which the growth of perturbations is suppressed are determined 
not only by the mass of the spherical component, but also by the density distribution in 
it. The suppressing effect is enhanced with increasing concentration toward the center. 

(2) The presence of a compact and moderately massive bulge in a galaxy effectively 
prevents the growth of bending perturbations. 

(3) It follows from an analysis of the entire set of our results that a more accurate 
approach to estimating the dark halo mass from the observed relative thickness of the 
stellar disk zo/h in spiral galaxies is required. 
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Fig. 1: (a) Magnitude of the vertical gravitational force (— -Fz) versus z at various distances 
R from the disk center; (b) the square of the vertical oscillation frequency v\ = —dF^/ dz 
versus z for various R. We took G — 1 and the disk parameters = 1 (the total disk 
mass), h — 3.5, and Zq — 0.1. 
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Fig. 2: Ratio <Jz/(Jr versus R for the time t = 3000 (model 12 with a massive halo: 
H = 3.0, lib = 0). 
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Fig. 3: Initial rotation curves for models 50, 76, 75, 49, and 53. The ratio of the total 
mass of the spherical component to the mass of the disk within a radius of 4h is the same 
for all models, = 0.5; /Xb = for model 50, /Xb = 0.0625 for model 76, /Xb = 0.125 for 
model 75, //b = 0.25 for model 49, and //b = 0.5 for model 53. The unit of velocity is 
978 km s-^ 
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Fig. 4: Evolution of the velocity dispersions and cr^ at i? = 2h for the same models 
those in Fig. El 
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Fig. 5: Edge-on view of the galaxy at the time t = 3000 for models 50, 76, 75, 49, and 53. 
The strenght of the image blackening corresponds to the logarithm of the particle number 
per pixel. The horizontal and vertical scales are 60 and 10, respectively. The ratio of the 
total mass of the spherical component to the mass of the disk within a radius of Ah is the 
same for all models, ii — 0.5. 
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Fig. 6: Radial thickness profiles for the galaxy at the time t — 3000 for model 53. The 
thickness was determined by two methods: as the rms value of the z coordinates of the 
disk particles, Zrms, ^-nd as twice the median of \z\, 2zi/2- 
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Fig. 8: Evolution of the disk thickness (2^1/2) at i? = 2/i for models 50, 76, 75, 49, and 53 
(see the caption to Fig. El). 
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Fig. 9: Radial thickness profiles for the galaxy (22:1/2) at the time t = 3000 for models 
that differ only by the initial thickness. For all of the models, jj, ~ 0.6 and /i-b = 0. 
Models 11_1, 11-2, and 11_3 are different random reahzations of a stellar system with 
Zq = 0.1; models 25 and 25_1 are different random realizations of a system with zq — 0.2; 
model 26_1 is for zq — 0.3. 
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Fig. 10: Radial thickness profile for the galaxy (22:1/2) at the time t = 3000 for models 
unstable against the growth of a bar mode {Qt — 1-5): (a) models with /j, — 0.5 (model 52 
with /^b = 0, model 51 with //b = 0.25), (b) models with fj, = 0.875 (model 43 with 
Hh = 0.25, model 42 with //b = 0.25); the initial half-thickness for all of the models is 
zo = 0.1. 
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